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ABSTRACT 

We introduce a Monte-Carlo simulation approach to thermodynamic Bethe 
ansatz (TBA). We exemplify the method on one particle integrable models, which 
include a free boson and a free fermions systems along with the scaling Lee- Yang 
model (SLYM). It is confirmed that the central charges and energies are correct to 
a very good precision, typically 0.1% or so. The advantage of the method is that 
it enables the calculation of all the dimensions and even the particular partition 
function. 



Conformal field theory in two dimensions and its perturbed integrable models 
have attracted considerable attention, since the work [1]. A great deal of interest 
stems from the work of Al. Zamolodchikov who considered first the thermalization 
of integrable systems (see ref. [2] and ref. therein). See also [3], and for further 
developments see [4]. Here we offer an approach to thermodynamic Bethe ansatz 
based on a Monte-Carlo simulation. 

We consider a a collection of particles moving in a box of length /. For simplicity 
we assume only one particle specie in a purely inelastic scattering matrix, S, and 
also assume periodic boundary condition on the box. 

We assume that the particles are at the locations xi,X2, ■ ■ ■ x n and that the 
distances X{ — Xj are large, so that we can use the scattering matrix to compute 
the relative phase. We denote the scattering amplitude by S(f3) where (3 is the 
relative rapidity, p = m sinh (3 is the momentum, E = m cosh (3 is the energy, and 
where m is the mass of the particle. 

The Bethe ansatz assumes that the wave function is 



where the theta ensures that x\ < x 2 < . . . < x n and the dots correspond to other 
arrangements of the coordinates, picking up a factor of S(f3i — (3j) every time X{ 
crosses Xj. For example for two particles, the wave function is 



ip(xi,x 2 , ...,x n 



) = Y[ e iXrPr 9{x 1 <x 2 < ■■■ <x n ) + 



(1) 



r 



4>(xi,x 2 ) = e ix ^ +ix ^ [6(x2 - Xl ) + S{(3i - f3 2 )6{x 1 - x 2 )\ . 



(2) 



The amplitude in a unitary theory can be written as, 



S(0) = e iS W, 



(3) 



where 8{(3) is a real phase. 



2 



The basic equation follows from taking a particle all around the box, picking 
up a phase each time we exchange two particles location and a phase from the 
wave function factor e iPiXi . It is 

3 
3^i 

Taking a log of this equation we arrive at the Bethe equation, 

ml sinh A + ^<KA - Pj) = 27m;, (5) 

3 

where i and j take the values 1 to n, and n is the number of particles, n; are some 
integers which determine the energy levels* 

The energy of this collection of particles is given by 



E = rE = mr cosh (6) 
i=l 

where r is the inverse temperature, 

where T is the temperature, and k is Boltzmann constant. We redefined the energy 
by multiplying it by r and denoted by E. 



* Denoting the wave function for these rij's by tp ni ,n 2 ,...,n m (xi), the actual wave function has 
to take into account that these are identical particles. Thus, the full wave function assumes 
the form: 

i>( x i) = Xl^ (p) ^p(«i),p(«2),...,p(« m )(^)' 

V 

where p is any permutation of m objects, and s(p) is (1) if p is an even (odd) permutation, 
respectively, and h = 1 (h = —1) for Bosc (Fermi) statistics, respectively. It follows that the 
wave function is defined for n\ < ri2 < . . . < n m for Bose statistics, and, ni < ni < . . . < n m 
for Fermi statistics. 
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We define a grand canonical partition function by summing over all the particle 
numbers, m, without a chemical potential, \x = 0, and summing over the energy 
levels, rii, taking into account that all the particles are identical. These means that 
we have to limit the sum on n\ < 112 < . . . < n m . The partition function then 
becomes, 

Z(m,r,l)= exp[-^("i)]- (8) 

,m 

n 1 <n 2 <...<n m 

In the case of Fermi statistics, when the ni are all diferent we may write the 
partition function as, 

Z(m,r,l) = J2 ^exp[-E(rn)]. (9) 

,m 

The average energy is given, as usual, by 

d 

E(m,r,l) = (E) = -r— log Z(m, r, /). (10) 

or 

The conformal limit is obtained by taking the mass to zero, m — > 0. In this 
limit the partition function becomes only a function of the ratio r = r/l. 

lim Z(m,r,l) = Z(r/l), (11) 

m— >0 

and we denote by E{r/l) = lim m ^o E(m, r, I) the average energy in this limit. We 
expect the partition function to be equal to the path integral on the torus with a 
modulus it of the corresponding conformal field theory. Denote by H the Hilbert 
space of this theory. Then the following correspondence emerges, 

Z(r/l) =Tre- 27rr ( Lo+Z °)/', (12) 

TC 

where Lq and Lq are the left and right moments of the stress energy tensor, i.e., 
the Hamiltonians of the system, whose eigenvalues are the conformal dimensions. 
The relation, eq. (12), can be used to evaluate all the dimension of the CFT and 
even the particular modular invariant used, by evaluating the partition function. 
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We can now use modular invariance [5, 6]. This means that the function, 

f{r/l) = Z(r/l)e 2nrc ^ 12l \ (13) 

where c is the central charge, is invariant under / and r replaced, 

f(l/r) = f(r/l). (14) 

By deriving the log of this equation with respect to r, we get a relation among 
energies, 

E(r/l) + E(l/r) = 2 -^-(L + ! j y (15) 

This is a very useful relation for evaluating the central charge. An exception to 
this relation occurs for free bosons where: 1) There is an infinite factor in Z due 
to the zero mode, ni = 0, which we can ignore (or eliminate by giving a chemical 
potential /z, for the rii = mode. 2) / has to be multiplied by 

(l/r^fir/l), (16) 

to make it modular invariant. This changes the relation eq. (15), to be 

E(r/l) + E(l/r) + 1 = ^ (I + ^ , (17) 

where the central charge c = 1 for one free boson. 

We start by studying some simple examples, which we take to be systems of 
one free boson or one free fermion. 
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For a free boson the different modes are independent and the partition function 
for the system is readily evaluated to be, 

oo oo oo 

Z(r/l) = U e~ 2nr ^ Sm / 1 = - q n y 2 , (18) 

s m =0 n=l 

where s m is the number of particles at the mth energy level and we defined q to be 

q = e- 2nr / 1 . (19) 

We ignored the zero mode which is a constant (infinite) factor, which can be 
controlled by adding a chemical potential, as explained above. 

We see immediately, that eq. (12) is indeed obeyed, and that Z is exactly the 
partition function of one free boson, on a torus with a modulus it. We can also 
verify that equation (17) holds giving the central charge c = 1 as expected. 

Let us consider now a system of one free fermion. We assume that 5(f3) = 
and that the boundary condition is anti-periodic. We take a Fermi statistics, i.e., 
only one particle is allowed for each energy level. This implies that the Bethe 
equation, eq. (5), becomes, 

mlsmhfii = 2vr(n; + 1/2), (20) 
and the partition function becomes, 

oo 

Z(r/l)= l[(l + q m+ h 2 , (21) 

where q is again given by eq. (19). 

Again, this is precisely the partition function for a free fermion in the Neveu- 
Schwarz sector, confirming eq. (12). We can verify that eq. (15) holds, giving the 
correct central charge, c = 1/2. 



6 



The Ramond sector is obtained by taking the trivial scattering matrix S — 1, 
periodic boundary conditions and Fermi statistics. The partition function can be 
easily evaluated using eqs. (5,9), and it is 



Z{r/l) = 2\[{l + q n ) 2 . (22) 



n=l 



We see that, indeed, it is identical to the free fermion CFT partition function, 
which is 



Z ctt (r/l) = 2q 1 ' 12 l[(l+q n ) 2 - (23) 

71=1 

We note however, that the overall factor is different. Thus we define 

XR = (1~ A+C/M XR, (24) 

where A is the dimension of the primary field corresponding to this block, and 
Xr is the usual Ramond character. We then see that the correspondence with the 
CFT is, 

Z(r/l) = 2\ XR \ 2 . (25) 



We conclude that two differences arise, in general, in eq. (12): 1) Each block 
in the CFT occurs for a different boundary conditions. 2) The overall factor eq. 
(24) has to be eliminated, making the character into a "q-series", i.e., only integer 
non-negative powers of q are allowed in the character. 

Now, for an interacting system it is quite difficult, in general, to compute 
this partition function. For example, the Lee- Yang theory M(2,5), which is the 
perturbed minimal model with p = 2 and q — 5 [1]. This theory has only one 
particle in the spectrum, particle A with a mass m and the scattering amplitude 
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is [7], 



/m sinh/3 + isin(7r/3) 
smh p — i sm(7r/3) 



Thus, we resort to a numerical algorithm for evaluating the partition function 
Z(m,r,l). There are two stages in this algorithm, which are basically, 1) Solving 
the basic equation (5). 2) Simulating the grand canonical ensemble eq. (8). 

To solve the equation we must calculate E(ni) for any set of levels, n{. This 
we do by taking the initial guess (3f = and solving eq. (5) by iterations, 



Pf = sinh" 1 



ml 



(27) 



where s denotes the sth iteration. Usually 5 or so iterations are enough to get 
with sufficient precision (10^ 6 or so). 

The second stage we use is to make a Monte-Carlo simulation of the system 
using a Metropolis type algorithm, in order to find the average energy, E(m,r,l). 
We use three types of steps, preserving the detailed balance to get the correct 
population of particles. 

For the first step, we change any of the n s 's for all the s by one n s — > n s ± 1, 

ni, n 2 , . • . , n m — ► ni, n 2 , . . . , n s -\,n s ± 1, n s+1 , . . . , n m . 

If the energy of the new configuration is less, we accept the new configuration. If 
the energy is greater, we accept it with the probability 

exp[i? i d - -Bnew], (28) 

which ensures the detailed balance distribution of probability exp[— E\. This step 
allows us to move in momentum space, corresponding to a thermalized random 
walk in the n lattice, covering all the possible values of n s (ergodicity). 
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The second and third steps allow for changing the numbers of particles, i.e., 
we either insert a new level somewhere, n s = for some s, or eliminating a level 
which is zero n r = 0. The adding of a particle corresponds to the step: 

ni, n 2 , ■ • . , n m — > ni, n 2 , . . . , n s _i, 0, n s+ i, . . . , n m (29) 

and the elimination of a particle corresponds to the step 

ni, n2, . • • , n s -i, 0, n s+ i, . . . , n m — > m, n 2 , ■ ■ ■ , n m . (30) 

Again we take the new configuration if its energy is lower. If the new con- 
figuration is higher in energy, we take the new configuration with the probability 
eq. (28). This allows the change of the number of particles and the population of 
energy levels is according to the grand canonical distribution eq. (8). The relation 
eq. (28) ensures the correct detailed balance. We use configurations of n s = as 
"gateways" to change the number of particles, deleting or adding only zero energy 
levels. 

Actually, the algorithm above does not take into account the ordering of the 
energy levels, n\ < 112 < . . . < n m . This we do in a different way for bosons and 
fermions. For fermions we add a particle with the probability exp[— AE]/(m + 
1) if the energy is higher, and 1 if it is lower. We delete a particle with the 
energy exp[— AE] if the energy is higher and l/(m + 1) if it is lower. If the ni 
already contains some n-i = we do not add a particle, and leave the configuration 
unchanged. Likewise if the configuration does not contain we do not delete a 
particle. 

For free fermions this algorithm can be simplified since adding a particle always 
raises the energy, and deleting it always lowers it. So, we add a particle with the 
probability exp[— AE]/(m + 1) and remove it with the probability 1. 

For bosons we use the same algorithm but allow only the configurations which 
are ordered: n\ < n 2 < . . . < n m discarding any other new configurations. We 
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change rii — > rn ± 1 at some location. If the new configuration is not ordered, we 
do not make the change. This guarantees detailed balance. 

This allows us to calculate the energies, eq. (10), preserving correctly the 
detailed balance. 

The easiest systems to simulate are, of course, the free fermions and free boson 
systems. For a free fermion system we took 8{f3) = along with antiperiodic 
boundary condition and Fermi statistics. For r/l < 1 we took 10000000 sweeps 
for each temperature, and for r/l > 1 we took 400000000 sweeps. The results are 
listed in tables (1) and (2). Here E mc and iV mc are the average energy and average 
particle numbers, as found in the simulation. We define the "partial" central charge 
as: 

c (r/l)- 12 ^™fr/0 (31) 
Cmc[r/l) ~ 2n(l/r + r/lY W 

and the central charge is given by, using eq. (15), 

Cmcir/l) + c rnc (//r) = c, (32) 

where c is the Monte-Carlo value of the central charge of the theory. For example 
we can calculate, from tables (1-2), 

c mc (0.5) + c mc (2) = 0.481836 + 0.0178107 = 0.499647, (33) 

which is impressively close to the actual central charge c = 0.5, in accordance with 
eq. (15). 

E'calc is the energy calculated from eq. (10). Here the average number of 
particles iV ca i c is calculated by introducing a chemical potential /i, 

Z(r/l,n)= exp[-E(rn) ~ »m], (34) 

rij ,m 
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for Bose statistics. For Fermi statistics we have, 



Z(r/l,/i)— exp[— E{rii) — jj,m] = — - exp[— E{rii) — fj,m]. (35) 

For a free fermion Z assumes the form, 

oo 

Z(r/l, fi) = ] J (1 + exp[-27rr(n + 1/2)// - /i]) 2 , (36) 
and for a free boson it is 

00 

Z(r/l, fi) = Y[ (1 - exp[-27rrn// - fi])~ 2 . (37) 
n=l 

The energy is given as in eq. (10), 

£ C alc = -r^logZ(r/Z,0), (38) 
We calculate the derivative with respect to ji to get the average number of particles, 

iV calc = - lim A log Z(r/l, /i). (39) 

H^O O/jL 

From tables (1-2), we see that across the range of temperatures the average 
energy and average iV are typically only 0.1% off the calculated numbers, which is 
a very reasonable correspondence. This demonstrates the very agreeable efficacy 
of the Monte-Carlo approach. 

In tables (3-4) we list the results for a free boson system, for r// < 1 (table 
(3)) and r/l > 1 (table (4)). We take a chemical potential for the zero modes, 
fi = 0.5. At this value only a few rii = are created. Also we make 100000000 
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sweeps of the Monte-Carlo. Again, the correspondence between the Monte-Carlo 
results and the calculated ones is very good, around 0.1%. We defined here c mc as 

12[S(r/Q + l/2] 

c mc {r/i)- 2n{l/r + r/l) ■ m 

It follows from eq. (17) that the Monte Carlo value of the central charge is, again, 

c mc (r/l) + c mc (//r) = c mc . (41) 
For example, from tables (3-4), 

c mc (0.8) + c mc (1.25) = 0.528855 + 0.471513 = 1.000368, (42) 

which is very close to the theoretical value c = 1. 

The zero modes of the bosonic system, i.e., for ni = have to be treated 
carefully, since they have zero energy. First, we introduce a chemical potential, 
only for the zero modes, in the Monte-Carlo simulation. Second, in calculating 
N mc we count only the number of the non-zero modes, in order to correspond with 
eq. (37), the partition function. We see from tables (3-4) that this leads to a very 
good agreement with the theoretical results. 

We now turn our attention to the scaling Lee Yang model, which is the per- 
turbed minimal model M(2,5). We first need to establish some facts about the 
minimal models, in general [1]. The models are labeled by two integers p and p' 
assumed to be strange to each other, and labeled by M(p,p'). The central charge 
is 

c=l-^2! (43) 
pp' 

The fields are given by n = 1, 2, . . . , p — 1 and m — 1, 2, . . . , p' — 1, denoted by n ,m 
whose dimensions are 

A {np 1 -mpf - {p-p'f ( 

ZVm ~ 4pp> ■ lMj 

The field <p 

n,m is identical to the field (J)p—n,p'—m- These models are unitary only for 
p' — p + 1 [8]. In particular, the model M(2, 5) is not unitary. However, it is only 
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"weakly" non-unitary, since the scattering matrix Saa, eq. (26) is unitary and 
also, as we shall see, the model is modular invariant. This enables us to calculate 
the partition function as if it was a unitary theory. 

The character of the field 4> m ^ n is defined by 

Xn ,m = Jr exp[2mr(L - c/24)}. (45) 
We define the classical theta functions at the level p by [9], 

©n, P (r)= £ (46) 

j=^-+integer 

where n is defined modulo 2p. The theta functions are modular forms, transforming 
by the modular transformation r ^ -1/r as [9] 

6 n , p (-l/r) = (-it)- 1 ' 2 £ ^e^(r), (47) 

n 

where the sum is over n modulo 2p and the matrix S is given by 

S n , n = J=e-" im /r. (48) 



The charaters of the minimal model, M(p,p') can then be seen to be, 

Xn,m(r) = en " Pp,(r) ?? ~ ) 9n+ ' Pp,(r) , (49) 

where we defined, 

n± = np' ± mp, (50) 
and where the Dedekind eta function is defined by, 

oo 

v(r)=q^ 2A ll(l-q n ). (51) 
n=l 

It is immediately seen that these characters give the correct central charge and 
dimensions, eq. (43-44). Moreover, they give the null vector for the field § n ,m at 
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the levels A = A niJn + mn and A = A n ^ m + {p — n){p' — m), in accordance with 
[!]■ 

The characters, so defined, are modular functions of a subgroup of the modular 
group. They are seen to transform according to, 

Xn,m(-l/r) = ^^ njm ,n,mXn,m(r), (52) 

where the matrix FT is can be seen, from eq. (48), to be 

HV*, - ^(-D— sin sin {^f) . (53) 

It is immediately seen that the theory M(p,p') is equivalent, at the level of the 
modular matrix, to SU(2) x SU(2) at the "pseudo" levels p/p' and p' jp [10, 11], 

M(p, p') » S£/(2) p/p , x 5C/(2 (54) 

The factor (_i) nm + jnn j n equation (53) can be ignored since the left movers and 
the right movers differ by an even integer. Here, twice the SU (2) isospin is given 
by j = ra — l, / — m — 1, j =7i — l, / =7ra — l, and the level is k = p — 2, k' — p' — 2. 
In this notation, the modular matrix of SU(2) p ,/ p is given by [9, 6, 10, 11] 



and, indeed, eq. (53) is exactly a product of two such factors, for p/p' and for p 1 } jp. 
For the modular transformation T, which is generated by r — > r + 1, we have 

T nm = e 2vri(A n , m -c/24) _ e 27rinV/(4p) e 27rim 2 p/( V) ^_]_^mn e -27ric/24 

The factor (— l) mn is irrelevant since it cancels between the right and left movers, 
as they differ by an even number, and so is the factor with the central charge. The 
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other factors correspond exactly to the dimensions of SU(2) p / p , and SU(2) p ,/ p , 

A D / j(j + 2) (57) 

and 

We infer that, also, the modular matrix T is a product of two S77(2)'s at the 
pseudo levels p and p' . (At the level of the dimensions, the correspondence with 
SU{2) x SU(2) was noted already in [12]). 

We conclude that the general modular invariant partition function is given by, 

Z ^ = \ N n,n'K mtm 'Xn,m(r)Xn',m'(ry, (59) 

n,m,n' ,m' 

where the factor of 1/2 accounts for the field identifications, and where N and 
K are any modular invariants of 577(2) at the levels k = p — 2 and k' — p' — 2, 
respectively, which are in relation with the simply laced Lie algebras of types ADE 
[6, 13]. For a proof see [14]. This solves the problem of classifying the acceptable 
partition functions. 

The modular transformations also imply that the fusion rules of the model 
M(p,p') are the same as a product of two S77(2)'s, 

4>n,m X <\)n,fh = ^ fn,n,qfm,fh,t ( i ) q,t-, (60) 
q,t 

where f^n, q is the fusion rule of SU (2) at the level k = p — 2 according to the 
"depth rule" [6] . We have that f% nq is equal to one if: 

n+n+g = 1 mod 2, q > \n—n\ + l, q < min(n+n— 1, 2p— n— n— 1), (61) 

and is zero otherwise. These fusion rules can be seen to be in accordance with the 
Verlinde formula [15]. When we use the S matrix eq. (53) we recover the correct 
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fusion rules, eqs. (60-61). This implies that these are fully consistent conformal 
data for a unitary theory, obeying all the axioms of unitary conformal data. 

Actually, the model M(p,p'), for p — p' odd, has exactly the same conformal 
data as the unitary model 

L = SU(2) 1 x SU(2) (pl/p) _ 1 x SU(2) {p/pl) _ 1: (62) 

where we include in the chiral algebra the field p — [1] x \p — 2] x [p' — 2], and 
where each of the numbers correspond to twice the isospin of SU(2). The central 
charge of the model is calculated to be, 

c = 1 + 3(p7p - l)(p - 2) + 3(p/p' - l)(p' - 2) = 1 - 6(P ~ p f )2 , (63) 

i.e., exactly as the central charge of M(p,p'), eq. (43). The field p has integer 
dimension, for p — p' odd, and so we can include it in the chiral algebra. Also, the 
consistency of the pseudo conformal field theory, i.e., at the level p — p' requires 
that p — p' is odd [10]. The extended field implies: 1) The field identifications, 
<pn,m = 4>p-n,p'-m, by multiplying any of the fields with this field. 2) The total 
isospin has to be an integer, for locality with respect to this field: [s] x [j] x [/] 
obeys s + j '• + I = mod 2. It can be easily verified that the modular properties 
of the theory L, with this extended field, are exactly the same as M(p,p') theory, 
eqs. (53,56). In particular, this implies the fusion rules, eqs. (60-61), using, for 
example, the Verlinde formula [15]. For p' = p + 1 we recover the conformal data 
of the coset model SU(2) 1 x SU(2) p / SU(2) p+1 , as is well known [16]. 

A similar correspondence holds when p and p' are both odd. In this case we 
define the model 

M » SU(2) (p+pl)/p /Z 2 x SU(2) p/p ,, (64) 

where SU(2)k/ 'Z 2 stands for the conformal field theory of integer isospin represen- 
tations of SU(2), see [10]. It can be immediately verified that the central charge 
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is correct, giving eq. (43), up to an integer, and that the modular transformations 
are identical with eqs. (53,56). We conclude that for all the minimal models, their 
conformal data can be realized by full fledged unitary conformal field theories. 

For the model M(2,5), the SLYM, we have two fields = 1, which is the 
unit operator, and (f> = 0i ; 2 which has the dimension Ai ; 2 = —1/5. According to 
eq. (59), only the diagonal modular invariant is allowed in this case. 

The fusion rules are easily read from eqs. (60-61) and are given by, 

2 = 1 + 0. (65) 

These are the same fusion rules as (£2)1 and SU(2)z/ Z2. The later is the con- 
formal theory consisting of the integer isospin representations of SU(2)s, see [10]. 
However, the modular matrix of M(2, 5) is at pseudo level 3 of SU{2)%/ Z2 and so 
is different from these cases (which are at the pseudo levels ±1). 

It is convenient to define the characters, 

WM=? A '" +c/2 Vn(r), (66) 

so that x contains only integer powers of q (a "q-series") starting at 1. The square 
of this object is what we actually compute in the Monte-Carlo simulation, and 
with which we want to make the comparison. We have approximately, using eq. 
(49), the character of the identity, 

1 — q — Q + Q + ■ ■ ■ 

(i- 9 )(i -,»)(! -,.)...• m 



and of the field 0, 



1 — (f 1 — q^ + + 

XlM = T\ \f\ 3^ ' ( 68 ) 

(1 -q)[l - q 2 )[l - q 6 ) . . . 



where we defined, as usual, 

q = exp(27nr), (69) 
and were, for the comparison, r = ir/l, following eq. (12). The square of one of 
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these two quatities, we expect to recover in the Monte-Carlo simulation, according 
to eq. (12). 

Now, in the conformal limit, m — > 0, we can, actually, solve exactly the Bethe 
ansatz for the SLYM system, eq. (5). Inspecting the scattering matrix in this case, 
5aa(/3)i we see from eq. (26) that 

Saa{P) = e*"W , (70) 

and that for (5 — 0, 8aa(P) — 71 an d that for (5 » 1, 8aa{0) — 0- This means 
that in eq. (5), since we can assume that if /5, — f3j is not zero, it is very big, so if 
rii ^ rij then the Bethe equation for i and j decouples. This implies that we need 
to solve the partition function only for n^s which are all the same. In this case, we 
simply need to take n{ — » rii — 1/2 if the number of particles is even, and to leave 
it unchanged, rii, if the number of particles is odd. 

We conclude that the partition function for SLYM is: 

oo 

Z(r/l) = H (1 + q s + q 2s ~ l + g 3s + q 4s ~ 2 + g 5s + g 6s - 3 + . . .f . (71) 
s=l 

(We omited the zero mode, which gives an infinite factor, and as explained before, 
it can be controlled by adding a chemical potential for the zero mode.) It is 
immediately observed that this partition function, eq. (71), is different from the 
square of the characters, eqs. (67-68). We conclude that the scattering matrix, eq. 
(26) needs to be modified by a "Z-factor" to make it correct for the SLYM system. 
We conjecture that it needs to be multiplied by the solution of the Yang Baxter 
equation based on the Hard Hexagon model, IRF{{G2)\,4>,4>) [10]- Further work 
on this is required. 

We hope that the Monte-Carlo simulation approach to the thermodynamic 
Bethe ansatz can be a valuable tool in the investigation of integrable systems, and 
as demonstrated here it can lead to very precise results. 
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Table 1. 
Free Fermion Results, r/l < 1. 



r/l 


E 


N mc 




-^calc 


A'calc 


0.05 


5.23606 


4.39815 


0.49876 


5.24908 


4.40616 


0.1 


2.64036 


2.19051 


0.49928 


2.64417 


2.19318 


0.15 


1.78325 


1.45019 


0.499623 


1.7846 


1.451 


0.2 


1.36164 


1.0766 


0.500103 


1.36135 


1.07633 


0.25 


1.11179 


0.847974 


0.499617 


1.11256 


0.848415 


0.3 


0.951349 


0.694035 


0.500076 


0.950611 


0.693539 


0.35 


0.838011 


0.580287 


0.499037 


0.837359 


0.580059 


0.4 


0.753898 


0.492667 


0.496496 


0.753123 


0.492233 


0.45 


0.686255 


0.421321 


0.490472 


0.686626 


0.421479 


0.5 


0.630722 


0.362648 


0.481836 


0.631075 


0.362824 


0.55 


0.583235 


0.313676 


0.470359 


0.582343 


0.313237 


0.6 


0.536727 


0.270202 


0.452238 


0.537975 


0.270783 


0.65 


0.496334 


0.234034 


0.433148 


0.496578 


0.234155 


0.7 


0.458426 


0.202888 


0.411322 


0.457424 


0.202418 


0.75 


0.419849 


0.174717 


0.384889 


0.420183 


0.174866 


0.8 


0.385041 


0.151055 


0.358719 


0.384752 


0.150935 


0.85 


0.351453 


0.13026 


0.331229 


0.351144 


0.130158 


0.9 


0.319961 


0.112321 


0.303853 


0.319412 


0.112135 


0.95 


0.289654 


0.0965256 


0.276236 


0.289617 


0.0965201 


I. 


0.262427 


0.0832056 


0.250599 


0.261799 


0.0830094 



Acknowledgements: I thank A. Babichenko, A. Belavin and S. Pal for helpful dis- 
cussions. 
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Table 2. 
Free Fermion Results, r/l > 1. 



r/l 




N mC 




-^calc 


-^calc 


1. 


0.262017 


0.0830768 


0.250208 


0.261799 


0.0830094 


1.05263 


0.234377 


0.0706784 


0.22352 


0.234671 


0.0707662 


1.11111 


0.206624 


0.0590791 


0.196221 


0.207096 


0.0592151 


1.17647 


0.179296 


0.0484496 


0.168978 


0.179385 


0.0484738 


1.25 


0.151874 


0.0386435 


0.141491 


0.151936 


0.0386596 


1.33333 


0.124918 


0.0298078 


0.114516 


0.125233 


0.0298832 


1.42857 


0.100163 


0.0223123 


0.089871 


0.0998352 


0.0222393 


1.53846 


0.0764025 


0.0158057 


0.0666761 


0.0763599 


0.015797 


1.66667 


0.0555191 


0.0106026 


0.0467796 


0.0554371 


0.0105871 


1.81818 


0.0377313 


0.00660547 


0.0304291 


0.0376456 


0.00659048 


2. 


0.0233141 


0.0037105 


0.0178107 


0.0234235 


0.00372794 
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Table 3. 
Free Boson Results, r/l < 1. 



r/l 








-^calc 


-^calc 


0.1 


4.29141 


3.80843 


0.906032 


4.28835 


3.80781 


0.15 


2.56987 


1.83794 


0.860102 


2.5692 


1.8375 


0.2 


1.72245 


1.03771 


0.816261 


1.72271 


1.03759 


0.25 


1.22361 


0.637259 


0.774554 


1.22529 


0.638052 


0.3 


0.903667 


0.414006 


0.737837 


0.902409 


0.41351 


0.35 


0.67905 


0.277354 


0.702127 


0.679256 


0.277446 


0.4 


0.517753 


0.190424 


0.670264 


0.518432 


0.190651 


0.45 


0.399507 


0.1334 


0.642885 


0.399148 


0.133236 


0.5 


0.308222 


0.0940281 


0.617436 


0.308909 


0.094242 


0.55 


0.239491 


0.0671612 


0.596375 


0.239727 


0.067244 


0.6 


0.186348 


0.0483229 


0.578307 


0.186231 


0.0482849 


0.65 


0.144824 


0.0348669 


0.562735 


0.14465 


0.0348313 


0.7 


0.112077 


0.025171 


0.549186 


0.112247 


0.0252108 


0.75 


0.0866057 


0.0182176 


0.537761 


0.086975 


0.0182923 


0.8 


0.0676613 


0.0133725 


0.528855 


0.0672725 


0.0132962 


0.85 


0.0519823 


0.00968116 


0.520219 


0.0519312 


0.00967727 


0.9 


0.0398549 


0.00702419 


0.512675 


0.0400067 


0.00705006 


0.95 


0.0308441 


0.00514874 


0.506253 


0.0307573 


0.00513968 


1. 


0.0235342 


0.00373636 


0.499938 


0.0235988 


0.00374886 
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Table 4. 
Free Boson Results, r/l > 1. 



r/l 




N m c 




-^calc 


-^calc 


1. 


0.0234887 


o 


00373087 

\J \J 1 \J W 1 


0.499895 


0.0235988 


0.00374886 


1.05263 


0.0177955 





00268717 


0.493809 


0.0178181 


0.00269044 


1.11111 


0.0130254 





00186362 


0.487197 


0.0130086 


0.00186161 


1.17647 


0.00916301 





00123849 


0.479864 


0.00912634 


0.00123387 


1.25 


0.00611186 





00077747 


0.471513 


0.00610499 


0.00077701 


1.33333 


0.00386213 





00046036 


0.461907 


0.00385577 


0.000460143 


1.42857 


0.00225915 





00025103 


0.450652 


0.00227017 


0.000252884 


1.53846 


0.00125341 





00012921 


0.437441 


0.00122546 


0.000126766 


1.66667 


0.000561796 





00005319 


0.421766 


0.000593164 


0.0000566413 


1.81818 


0.000261345 





00002242 


0.403444 
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2. 


0.000104389 
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0.00000697 
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APPENDIX 

Calculating the average energy and number of particles. 

Enclosed is a Mathematica program which calculates the partition function 
ff[a,l,mu], for one free fermion in a box of length /, the inverse temperature a 
and a chemical potential mu, using eq. (21). Then it calculates the average energy 
ee[a, I] and the average number of particles, nn[a,l}. For free bosons we use the 
function ffb[a,l,mu] listed below too, which gives the partition function of one 
free boson, without the zero mode, and with a chemical potential mu. 

ff[a_ ,L ,mu_ ]:=Product[(l+Exp[-mu-2 Pi a/1 (n + 1/2)]) A 2,{n,0,60}] 

ee[a_ , 1 _] := -x D[Log[ff[x,l,0.]],x] /. x - > a 

nn[a_ ,L ]:=-D[Log[ff[a,l,mu]],mu] /. mu — > 0. 

Print[ee[0.05,l]," " ,nn[0.05,l]] 

5.24908 4.40616 

ffb[a_,L,mu_]:=Product[(l-Exp[-mu-2 Pi a/1 n])"2,{n,l,60}] 
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